Answering Questions; Gather Data, 2021-05

# Check requisite packages are installed.
packages <- c(
  "plotly", 
  "dplyr",
  "RMTRCode2"
)
for (pkg in packages) {
  library(pkg, character.only = TRUE)
}

package 㤼㸱plotly㤼㸲 was built under R version 4.0.5package 㤼㸱ggplot2㤼㸲 was built under R version 4.0.3package 㤼㸱dplyr㤼㸲 was built under R version 4.0.4

# Reserved Names
candidateData <- NULL
islandInteractionsOneEmptyTwoWhich <- NULL
islandInteractionsOneTwo <- NULL
islandInteractionsOneTwoWhich <- NULL
mats <- NULL
paramFrame <- NULL
plotScalingData <- NULL
pools <- NULL

Disentangling Effects on the Viking Data

Load Data

ellipsisApply <- function(..., FUN) {
  lapply(as.list(...), FUN)
}

load("LM1996-NumPoolCom-QDat-2021-05.RData")
# Stop if not all are not null
stopifnot(all(unlist(ellipsisApply(
  FUN = function(bool) {!is.null(bool)},
  candidateData, 
  islandInteractionsOneEmptyTwo,
  islandInteractionsOneEmptyTwoWhich,
  islandInteractionsOneTwo,
  islandInteractionsOneTwoWhich,
  mats,
  paramFrame,
  plotScalingData,
  pools
))))

plotScaling <- plotly::plot_ly(
  plotScalingData,
  x = ~Basals,
  y = ~Consumers,
  z = ~CommunitySize,
  color = ~Dataset,
  colors = c("red", "blue", "black")
)

plotScaling <- plotly::add_markers(plotScaling)

plotScaling <- plotly::layout(
  plotScaling,
  scene = list(
    xaxis = list(type = "log"),
    yaxis = list(type = "log"),
    camera = list(
      eye = list(
        x = -1.25, y = -1.25, z = .05
      )
    )
  )
)

plotScaling

# Check that the Two island and Three island scenarios are set-up the same.
stopifnot(unlist(lapply(islandInteractionsOneTwo, length)) == 
            unlist(lapply(islandInteractionsOneEmptyTwo, length)))
stopifnot(names(islandInteractionsOneTwo) == 
            names(islandInteractionsOneEmptyTwo))
# Check that the Which versions correspond correctly.
stopifnot(
  unlist(lapply(islandInteractionsOneTwoWhich, function(x) {
    length(RMTRCode2::CsvRowSplit(x))
  })) 
  == unlist(lapply(islandInteractionsOneTwo, function(x) {
    # We're like onions; we have LAYERS!
    lapply(x, function(y) {
      lapply(y, function(z) {
        sum(z > 1E-6) # How many "large" entries are there?
      })})}))
)
stopifnot(
  unlist(lapply(islandInteractionsOneEmptyTwoWhich, function(x) {
    length(RMTRCode2::CsvRowSplit(x))
  })) 
  == unlist(lapply(islandInteractionsOneEmptyTwo, function(x) {
    # We're like onions; we have LAYERS!
    lapply(x, function(y) {
      lapply(y, function(z) {
        sum(z > 1E-6) # How many "large" entries are there?
      })})}))
)

# Hybrids
# Create a count of how many times each entry will be repeated.
communitiesAllRepeater <- 5 * unlist(lapply(islandInteractionsOneTwo, length))

# Create template.
communitiesAll <- data.frame(
  CombnNum = rep(0, sum(communitiesAllRepeater)), # Should repeat all rows.
  Basals = 0,
  Consumers = 0,
  Dataset = "",
  DatasetID = 0,
  Communities = "",
  CommunitySize = 0,
  OtherSteadyStates = 0, # To be recalculated
  CommunityAbund = "",
  CommunityProd = 0,
  TotalID = "",
  # Additional Column!, 1 for direct assembly, 0 unused.
  IslandsUsed = rep(c(2,2,3,3,3), sum(communitiesAllRepeater)/5)
)

# Retrieve the rows used to make hybrids
communitiesAllProspects <- candidateData %>% dplyr::group_by(
  CombnNum, Basals, Consumers, Dataset, DatasetID, TotalID
) %>% dplyr::select(
  CombnNum:DatasetID, TotalID
) %>% dplyr::summarise(
  Count = dplyr::n(), .groups = "keep"
) %>% dplyr::filter(
  Count > 1
) %>% dplyr::select(
  -Count
) %>% dplyr::arrange(
  DatasetID, CombnNum
)

# Make sure that the names match.
stopifnot(communitiesAllProspects$TotalID == names(communitiesAllRepeater))

# Insert repetitions.
communitiesAll$CombnNum <- rep(communitiesAllProspects$CombnNum, communitiesAllRepeater)
communitiesAll$Basals <- rep(communitiesAllProspects$Basals, communitiesAllRepeater)
communitiesAll$Consumers <- rep(communitiesAllProspects$Consumers, communitiesAllRepeater)
communitiesAll$Dataset <- rep(communitiesAllProspects$Dataset, communitiesAllRepeater)
communitiesAll$DatasetID <- rep(communitiesAllProspects$DatasetID, communitiesAllRepeater)
communitiesAll$TotalID <- rep(communitiesAllProspects$TotalID, communitiesAllRepeater)

# To move over from the data.
# Communities = "",
# CommunityAbund = ""

communitiesAll[communitiesAll$IslandsUsed == 2, ]$Communities <- 
  islandInteractionsOneTwoWhich
communitiesAll[communitiesAll$IslandsUsed == 3, ]$Communities <- 
  islandInteractionsOneEmptyTwoWhich

communitiesAll[communitiesAll$IslandsUsed == 2, ]$CommunityAbund <- 
  # We're like onions; we have LAYERS!
  unlist(lapply(islandInteractionsOneTwo, function(x) {
    lapply(x, function(y) {
      lapply(y, function(z) {
        toString(z[z > 1E-6])
      })
    })
  }))
  
communitiesAll[communitiesAll$IslandsUsed == 3, ]$CommunityAbund <- 
  unlist(lapply(islandInteractionsOneEmptyTwo, function(x) {
    lapply(x, function(y) {
      lapply(y, function(z) {
        toString(z[z > 1E-6])
      })
    })
  }))

# To calculate from the data.
# CommunitySize = 0, # To be calculated from Communities.
# OtherSteadyStates = 0, # To be recalculated last after filtering
# CommunityProd = 0, # To be recalculated after Abund stored.
communitiesAll$CommunitySize <- unlist(lapply(
  communitiesAll$Communities, function(x) {
    length(RMTRCode2::CsvRowSplit(x))
  })) 

for (r in 1:nrow(communitiesAll)) {
  communitiesAll$CommunityProd[r] <- with(
    communitiesAll[r, ], 
    RMTRCode2::Productivity(
      Pool = pools[[DatasetID]][[CombnNum]], 
      InteractionMatrix = mats[[DatasetID]][[CombnNum]], 
      Community = Communities, 
      Populations = CommunityAbund
    )
  )
}

# Original systems
communitiesAll <- rbind(
  candidateData %>% dplyr::select(
    -CommunityFreq, -CommunitySeq
  ) %>% dplyr::mutate(
    IslandsUsed = 1
  ), 
  communitiesAll
)

# Treating the Productivity like one might treat a hash,
# if two rows with the same properties are assigned the same hash, 
# we only keep one. 
# One decimal place might be excessive, 
# but we can reflect on that if results down the line are not interesting.
# For the record though, it appears that it is a decently good approach at 
# removing effectively numerical duplicates.
# Not bothering, sort of, with IslandsUsed, since many times a community is
# reproduced on varying numbers of islands.

# communitiesAll <- communitiesAll %>% dplyr::mutate(
#   tempProd = round(CommunityProd, 2)
# ) %>% dplyr::distinct(
#   CombnNum, Basals, Consumers, Dataset, DatasetID, TotalID,
#   Communities, CommunitySize, tempProd, IslandsUsed,
#   .keep_all = TRUE
# ) %>% dplyr::select(
#   -tempProd
# )

communitiesAll <- communitiesAll %>% dplyr::mutate(
  tempProd = round(CommunityProd, 2)
) %>% dplyr::group_by(
  CombnNum, Basals, Consumers, Dataset, DatasetID, TotalID,
  Communities, CommunitySize, tempProd,
) %>% dplyr::summarise(
  CommunityAbund = dplyr::first(CommunityAbund),
  CommunityProd = dplyr::first(CommunityProd),
  IslandsUsed = toString(unique(IslandsUsed)),
  .groups = "drop"
) %>% dplyr::select(
  -tempProd
) %>% dplyr::group_by(
  CombnNum, Basals, Consumers, Dataset, DatasetID, TotalID
) %>% dplyr::mutate(
  OtherSteadyStates = dplyr::n() - 1,
  Islands1 = grepl(pattern = "1", IslandsUsed, fixed = TRUE) # Will be useful
)

Persistence of Hybrid Communities

The idea is straightforward: after allowing interactions between islands, for islands that are not the same as one of the original communities, isolate the island and check to see what happens.

communitiesHybrids <- communitiesAll %>% dplyr::filter(
  !Islands1
) %>% dplyr::select(-Islands1)

communitiesHybrids$AfterSepAbund <- ""
communitiesHybrids$AfterSepCommunity <- ""
communitiesHybrids$AfterSepCommunitySize <- 0
communitiesHybrids$AfterSepProduction <- 0
for (r in 1:nrow(communitiesHybrids)) {
  temp <- with(
    communitiesHybrids[r, ],
    {    
      temp <- RMTRCode2::CsvRowSplit(Communities)
      RMTRCode2::LawMorton1996_NumIntegration(
        A = mats[[DatasetID]][[CombnNum]][temp, temp],
        R = pools[[DatasetID]][[CombnNum]]$ReproductionRate[temp],
        X = RMTRCode2::CsvRowSplit(CommunityAbund), 
        OuterTimeStepSize = 3E4,
        Tolerance = 1E-6
      ) # retrieve the abundance over time matrix
    }
  ) 
  
  temp <- temp[nrow(temp), -1] # choose last row, remove time column.
  
  communitiesHybrids$AfterSepCommunity[r] <- toString(
    RMTRCode2::CsvRowSplit(communitiesHybrids$Communities[r])[which(temp > 1E-6)]
  )
  
  temp <- temp[which(temp > 1E-6)] # remove microfoxes.
  
  communitiesHybrids$AfterSepAbund[r] <- toString(temp)
  communitiesHybrids$AfterSepCommunitySize[r] <- length(temp)
  
  communitiesHybrids$AfterSepProduction[r] <- with(
    communitiesHybrids[r, ], 
    RMTRCode2::Productivity(
      Pool = pools[[DatasetID]][[CombnNum]], 
      InteractionMatrix = mats[[DatasetID]][[CombnNum]], 
      Community = AfterSepCommunity, 
      Populations = AfterSepAbund
    )
  )
}

communitiesHybrids <- communitiesHybrids %>% dplyr::mutate(
  Persists = AfterSepCommunity == Communities,
  ProdChange = AfterSepProduction - CommunityProd
)

So after running the dynamics for 3E4 time units (i.e. 3x the length of time the dynamics in the numerical assembly runs in between assembly steps and 1.5x the length of the island dynamics), the communities that persist are 6, 7, 8, 9, 10. Examining the communities themselves, they are all the same community, albeit with different starting points.

communitiesHybrids[communitiesHybrids$Persists, ]

An obvious follow-up question: how many of the communities that collapse do so to communities that we have not already seen?

communitiesHybrids <- communitiesHybrids %>% dplyr::mutate(
  AfterSepCommunityAlreadyPresent = AfterSepCommunity %in% communitiesAll$Communities
)
sum(!communitiesHybrids$AfterSepCommunityAlreadyPresent)

[1] 12

Consolidating down to unique ending states we have the following.

communitiesHybrids[
  !communitiesHybrids$AfterSepCommunityAlreadyPresent, 
  ] %>% dplyr::distinct(CombnNum, AfterSepCommunity, .keep_all = TRUE)

We will add these new states to our catalogue of communities from the experiments. We also take the abundance after separation if the community persists to better reflect steady-state conditions.

communitiesAll <- rbind(
  communitiesAll %>% dplyr::filter(
    Islands1 == TRUE
  ) %>% dplyr::mutate(
    HybridCollapse = FALSE, Persists = TRUE
  ),
  communitiesHybrids %>% dplyr::mutate(
    CommunityAbund = ifelse(Persists, AfterSepAbund, CommunityAbund),
    Islands1 = FALSE, HybridCollapse = FALSE,
  ) %>% dplyr::select(
    -AfterSepAbund, -AfterSepCommunity, -AfterSepCommunitySize, 
    -AfterSepProduction, -ProdChange, -AfterSepCommunityAlreadyPresent
  ) ,
  with(
    communitiesHybrids[
      !communitiesHybrids$AfterSepCommunityAlreadyPresent, 
    ] %>% dplyr::distinct(CombnNum, AfterSepCommunity, .keep_all = TRUE),
    data.frame(
      CombnNum = CombnNum,
      Basals = Basals,
      Consumers = Consumers,
      Dataset = Dataset,
      DatasetID = DatasetID,
      TotalID = TotalID,
      Communities = AfterSepCommunity,
      CommunitySize = AfterSepCommunitySize,
      CommunityAbund = AfterSepAbund,
      CommunityProd = AfterSepProduction,
      IslandsUsed = IslandsUsed,
      OtherSteadyStates = 0,
      Islands1 = FALSE,
      HybridCollapse = TRUE,
      Persists = TRUE,
      stringsAsFactors = FALSE
    ))
)

Invadability of Hybrid Communities

Looking at a longer time scale, what happens if/when invasions resume? Do the hybrid communities that emerged retain the uninvadability of the parent communities?

This question should be straightforward as it is testing a step from the assembly process.

communitiesAll$Uninvadable <- NA
for (r in 1:nrow(communitiesAll)) {
  communitiesAll$Uninvadable[r] <- with(
    communitiesAll[r, ],
    {
      tempRow <- rep(NA, nrow(pools[[DatasetID]][[CombnNum]]) + 1)
      tempRow[RMTRCode2::CsvRowSplit(Communities) + 1] <- 
        RMTRCode2::CsvRowSplit(CommunityAbund)
      RMTRCode2::LawMorton1996_CheckUninvadable(
        AbundanceRow = tempRow,
        Pool = pools[[DatasetID]][[CombnNum]],
        CommunityMatrix = mats[[DatasetID]][[CombnNum]]
      )
    }
  )
}

We can compare this property against some of the other properties.

Uninvadability versus whether a community was found via assembly (“on Island 1”):

with(communitiesAll,
     table(Uninvadable, Islands1))

           Islands1
Uninvadable FALSE TRUE
      FALSE    20    0
      TRUE     32   30

Never invadable and assembled (good!), but about half of uninvadables are found without being assembled. What about of those that persist?

with(communitiesAll %>% dplyr::filter(Persists),
     table(Uninvadable, Islands1))

           Islands1
Uninvadable FALSE TRUE
      FALSE     8    0
      TRUE      0   30

Which of course fills in some of the blanks. So none of the communities that persist are uninvadable if they were not an end state of the assembly process.

Presence of Mass Effects

We check to see what happens when we treat each community as a pool for the other and perform assembly. Are the results the same as the diffusion system?

First, update the pairings.

communitiesAll <- communitiesAll %>% dplyr::group_by(
  CombnNum, Basals, Consumers, Dataset, DatasetID, TotalID
) %>% dplyr::mutate(
  OtherSteadyStates = dplyr::n() - 1
) %>% dplyr::ungroup()

This procedure can be done in two ways: first by directed invasion where one community is a pool for the other, and second with mutual (undirected) invasion where both communities are simultaneously pools for and invaded by each other. Note that in the directed case, we do not need to do any of the communities already marked as uninvadable with respect to the regional pools. The other communities they would be compared with are subsets of the regional pools, and so would already be checked against.

Indirect Mutualism

Here, we check to see if the networks created by each community (hybrid or otherwise) has mutualism embedded in it.